First principles investigation of finite— temperature behavior in small sodium clusters 
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A systematic and detailed investigation of the finite-temperature behavior of small sodium clus- 
ters, Nan, in the size range of n = 8 to 50 are carried out. The simulations are performed using 
density-functional molecular-dynamics with ultrasoft pseudopotentials. A number of thermody- 
namic indicators such as specific-heat, caloric curve, root-mean-square bond-length fiuctuation, 
deviation energy, etc. are calculated for each of the clusters. Size dependence of these indicators 
reveals several interesting features. The smallest clusters with n = 8 and 10, do not show any sig- 
nature of melting transition. With the increase in size, broad peak in the specific-heat is developed, 
which alternately for larger clusters evolves into a sharper one, indicating a solidlike to liquidlike 
transition. The melting temperatures show irregular pattern similar to experimentally observed one 
for larger clusters [ M. Schmidt et al., Nature (London) 393, 238 (1998) ]. The present calculations 
also reveal a remarkable size-sensitive effect in the size range of n = 40 to 55. While Na4o and Nass 
show well developed peaks in the specific-heat curve, Naso cluster exhibits a rather broad peak, 
indicating a poorly-defined melting transition. Such a feature has been experimentally observed for 
gallium and aluminum clusters [ G. A. Breaux et al, J. Am. Chem. Soc, 126, 8628 (2004); G. 
A.Breaux et al, Phys. Rev. Lett., 94, 173401 (2005) ]. 

PACS numbers: 61.46.-l-w, 36.40.-c, 36.40.Cg, 36.40.Ei 



I. INTRODUCTION 

Finite temperature studies of finite-sized systems have 
been a topic of considerable interest during last decades. 
Recent experimental as well as theoretical studies have 
brought out a number of intriguing results. In a series of 
experiments on free sodium clusters in the size range of 
55 to 350, Haberland and co-workers have observed a 
large size-dependent fiuctuation in the melting temper- 
atures. They have also observed a substantial lowering 
of about 30 % in the melting temperatures as compared 
to that of the bulk. Interestingly, recent experiments 
by Jarrold and coworkers 0, HI show that small clusters 
of Sn and Ga have higher-than-bulk melting tempera- 
tures. In our previous investigations, we have attributed 
this higher-than-bulk melting temperatures mainly to 
the covalent bonding in these clusters, as against metal- 
lic bonding in the bulk phase. 0,0 Very recently Breaux 
et al. have seen a remarkable size-sensitive feature of the 
melting transition of gallium as well as aluminum clus- 
ters. Pj'tI Their experiments show that the nature of the 
heat capacity curve changes dramatically with addition 
of few atoms. For instance, Ga.'^Q does not show obvi- 
ous melting transition, while Gag"]^ exhibit a well-defined 
peak, and Gaj^j shows a broad peak in the heat capacity. 

A number of computer simulations on melting of 
small sodium clusters have been reported in literature. 
Calvo and Spiegelmann Q have performed extensive 
simulations on Na clusters in size range of 8 to 147. 



•Electronic address; ImsieeQunipune. erne t . in| 
t Electronic address: chacko@ unipune . ern et . in] 
■f Electronic address: kanhereO unipune. ernet . in| 



Their simulations employed the second moment approx- 
imation (SMA) potential of Li et al. |^ as well as 
the distance-dependent tight-binding (DDTB or TB) 
method. They found more than one peak in the heat 
capacity most of of the clusters studied. They further 
observed that the nature of the ground-state geometry is 
crucial to precisely understand the thermodynamic prop- 
erties of clusters. However, although the method they 
employed provide relatively good statistics required to 
converge the features in the caloric curve, it did not in- 
corporate the essential ingredients of electronic structure 
effects. These simulations hence failed to reproduce the 
crucial features of the experimental results, clearly bring- 
ing out the importance of incorporating the electronic 
structures effects. In a very recent study 'lO^, we have 
successfully reproduced the melting temperatures of Na^r 
{N = 55,92, 142) using the Kohn-Sham (KS) based ap- 
proach 11] of the density-functional theory (DFT) and 
also gave a plausible explanation on its irregular varia- 
tion. There have also been a few theoretical investiga- 
tions on melting of sodium clusters with sizes n < 55. 
Manninen et al. ,12j have investigated the melting tran- 
sition of Na4o cluster using ab initio method. They raised 
the temperature of the system from 150 K to 400 K at 
the rate of 5 K/ps, and found that the melting transi- 
tion occurs at the temperatures between 300 ^ 350 K 
Na4o. They also show that Nag exhibits only isomeriza- 
tion. Aguado et al. 13] have performed Car-Parrinello 
orbital-free simulations to investigate melting phenom- 
ena in Nag and Na2o- Their simulation times were 8- 
60 ps per temperature. They observed a clear peak for 
Nag and double peaks for Na2o in the specific-heat curve. 
We have also investigated the melting transition in these 
clusters using various methods, and found a model de- 
pendence in the melting characteristics. [14| 
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In the present work, we perform density-functional 
molecular-dynamical simulations on Na„ clusters (n = 8, 
10, 13, 15, 20, 25, 40, and 50) to investigate their thermo- 
dynamic properties, specifically the size-dependent fea- 
tures. We perform simulations with about 150 ps per 
temperature which is much larger simulation times than 
any other earlier work. In addition to the standard indi- 
cators like specific-heat, caloric curve, root-mean-square 
bond-length fluctuations, we also calculate the energy 
deviation, the potential energy difference between the 
solidlike-state and the liquidlike-state, etc. 

In Sec. m we describe the computational details, fol- 
lowed by the results and discussion in Sec. IIIII Finally, 
we summarize the results in Sec. IIVI 

II. COMPUTATIONAL DETAILS 

We carry isokinetic Born-Oppenheimer molecular- 
dynamics calculations using Vanderbilt's ultrasoft 
pseudopotentials 'iG'l within the local-density approxi- 
mation (LDA), as implemented in the VASP package. 
We use two different methods to obtain the ground-state 
and several equilibrium geometries for each of the clus- 
ters. First, a "basin hopping" algorithm is employed 
to generate few tens of structures for smaller clusters 
and several hundreds structures for larger clusters us- 
ing the second moment approximation (SMA) parame- 
terized potential of Li et al. Q Several of these geome- 
tries, say the lowest 10-70 geometries, are then optimized 
using ab initio density-functional method. In the sec- 
ond method, we obtained few more equilibrium geome- 
tries by optimizing several structures selected from high- 
temperature ab initio molecular-dynamics runs, typically 
taken from temperatures near and well above the melt- 
ing temperatures of the clusters. The simulations have 
been carried out for 12 temperatures in the range of 
lOOK < T < 750K for n = 8 and 10, 9-12 tempera- 
tures in the range of lOOK < T < 4:50K for the rest 
clusters. For all the cases, the simulation time is 150 ps 
per temperature. We have discarded first 30 ps for each 
temperature to allow for thermalization. An energy cut- 
off of 3.6 Ry 20] is used for the plane wave expansion of 
the wavefunction, with a convergence in the total energy 
of the order of lO"** eV. The resulting ionic trajectory 
data have been used to study the melting of clusters by 
analyzing various thermodynamic indicators, which are 
discussed below in detail. 

We calculate the deformation parameter, e^e/, to an- 
alyze the shape of the ground-state geometry for all the 
clusters. The shape of the ground-state geometry plays 
a crucial role in determining the thermodynamic prop- 
erties of a cluster. The deformation parameter, Sdef, is 
defined as 

_ 2Qi 

''^^^ Q2 + Q3 

where Qi > Q2 ^ Qa are eigenvalues of the quadrupole 



tensor Qij = J^i RaRij with Rn being coordinate of 
ion / relative to the center of mass of the cluster. A spher- 
ical system {Qi = Q2 ~ Q3) has Sdef — 1, while Sdef > 
1 indicates a quadrupole deformation of some kind. 

To analyze the thermodynamic properties, we first cal- 
culate the ionic specific-heat and the average poten- 
tial energy per temperature (the caloric curve). We 
extracted the classical ionic density of states, i^{E), 
of the system, or equivalently the classical ionic en- 
tropy, S(E) = fcBlnf2(i?), via the multiple histogram 
method |22| to evaluate the canonical specific-heat. In 
the canonical ensemble, the specific-heat is defined as 
C(r) = dU{T)/dT, where UiT) = JE p{E,T)dE is the 
average total energy. The probability of observing an en- 
ergy E at a, temperature T is given by the Gibbs distri- 
bution p{E,T) = n{E)exp{~E/kBT)/Z{T), with Z{T) 
the normalizing canonical partition function. We nor- 
malize the calculated canonical specific-heat by the zero 
temperature classical limit of the rotational plus vibra- 
tional specific-heat, i.e. Cq = (37V — 9/2)fcs. Details of 
this method can be found in Ref. Q. The melting tem- 
perature has been taken as a peak in the specific-heat 
curve, following the convention of the experiments. 

Other thermodynamic indicator calculated is the root- 
mean-square bond length fluctuations (RMSBLF), i.e 
the Lindemann-like criterion for a finite system, given 
as 

2 ^ {R.)i 

where, Rij is the distance between i*'' and j*'* ion. This 
quantity gives the average fluctuation in the average bond 
lengths that are occurring at a given temperature. A 
value of about 0.1-0.15 signifies a melting transition for 
the bulk. However, as we shall see, for smaller clusters, 
this indicator should be taken with some caution. It is 
useful when examined in conjunction with other indica- 
tors such as the specific-heat. 

We have also calculated the energy deviation, 5E, |2^ 
defined as 

5E - {Etotai{T)) - [So + (3n - 6)fcsr] 

where, {Etotai{T)) = {Ek^T)) + {Epot{T)) is the aver- 
age total energy of the system at temperature T, with 
{Ekin{T)) and {Epot{T)) being the average kinetic en- 
ergy and the average potential energy, respectively. Eq 
is the ground-state energy and (3n — 6)kBT is the vibra- 
tional energy in the classical limit. 6E may be considered 
to be an indicator of anharmonicity in the system, as a 
function of temperature. Following Chuang et al. |23j| . 
we smoothened the plots of 6E using three point mov- 
ing average method. The error bars in the plot are the 
standard errors of every three data points. 

We have also examined carefully the role of simulation 
time. For this purpose, in Fig. ^ we plot the specific- 
heat for Na2o with two simulation times: one with 90 ps, 
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FIG. 1: The specific-heat of the Na2o cluster with simulation 
time of (a) 150 ps and (b) 90 ps as a function of temperature. 
Co = (3A'' — 9/2)fcs is the zero-temperature classical limit of 
the rotational plus vibrational canonical specific-heat. 

and another with 150 ps, per temperature. It may be im- 
mediately seen that the 90 ps data resuhs in a premelting 
feature which is absent for the 150 ps one. This indicates 
that even for smaller systems such as this, one needs to 
go to higher simulation times of the order of 150 ps or 
so. 



III. RESULTS AND DISCUSSION 

We begin by analyzing the geometries of sodium clus- 
ters for the sizes of n = 8, 10, 13, 15, 20, 25, 40 and 
50. This is then followed by a discussion on their ther- 
modynamic properties. We also address certain features 
of NaR5 and Na92 clusters relevant to the present discus- 
sion. [lO| 



A. Geometry 

The lowest-energy geometries of sodium clusters are 
shown in Fig. [3 First, we note that the equilibrium ge- 
ometries of Nas, Nai3 and Na2o, obtained by us, are in 
agreement with those reported by Rothlisberger et al. 
The ground-state geometry of Nag (Fig.|2Ia)) is a dodec- 
ahedron. One of its low-energy isomer is an antiprism 
(A E = 0.04 eV). These two structures play a crucial 
role in the finite temperature behavior of this cluster. We 
find two nearly degenerate structures for Naio, namely a 
bicapped dodecahedron (Fig.|2Ib-i)) and a bicapped an- 
tiprism (Fig.|2b-u)). Rothlisberger et al. |2^ have found 
the bicapped dodecahedron to be unstable. However, 
we computed the vibrational spectra for both geometries 
and found the structures to be stable. The clusters, Naia 
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FIG. 2: The ground-state geometries of the Na,i (n = 8 — 55) 
clusters. 
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FIG. 3: The deformation parameter, Edef, of the ground- 
state geometries as a function of cluster size. 

(Fig.[2Ic)) and Nais (Fig.j^^d)), exhibit capped pentago- 
nal bipyramidal structures as their lowest-energy config- 
urations. The ground-state geometry of Na2o (Fig.[2Je)) 
is a fivefold capped icosahedron, with two atoms capping 
the icosahedral faces on the central plane. We find a 
capped double icosahedron to be about 0.11 eV higher in 
energy than the ground-state. The lowest-energy struc- 
ture of Na25 has not been previously reported. It may be 
described as capped double icosahedron with growth on 



4 




8 10 13 15 20 25 40 50 55 

Size (No. of Atoms) 

FIG. 4: The eigenvalue spectra of their ground-state of Na„ 
(n — 8 — 50) cluster as a function cluster size. The numbers 
on the right show the degeneracy of that level. 



one side, leading to a distorted non-spherical structure, 
as shown in Fig. [2If). The ground-state geometry for 
Na4o, shown in Fig. [Hg), agrees with the one reported 
by Manninen et al. jl^] It consists of three decahedra 
capped by rest of the atoms. The structure is compact 
and retains the five-fold symmetry. The ground-state ge- 
ometry of Na5o, shown in Fig.l^fh), is highly asymmetric. 
It can be seen from figure [21 that a growth towards 50- 
atom structure starting from symmetric Na4o makes the 
structure non-spherical and asymmetric, which is simi- 
lar as seen for Na25 cluster. Finally, the ground-state 
geometry of Nass is a slightly distorted double Mackay 
icosahedron. It is the most spherical structure and is 
noted for the sake of completeness. ||lO| 

We have also examined the eigenvalue spectra and the 
shapes of the ground-state geometries of these clusters, 
as shown in Figs. 01 and |3| respectively. The shape defor- 
mation parameter, Edef, plotted in Fig.O for the ground- 
state geometries of all the clusters show that Na2o, Na4o 
and Na55 are nearly spherical, while Na25 and Naso are 
deformed. This behavior is also reflected in the eigen- 
value spectra of these clusters. The eigenvalue spectra of 
Na2o, Na4o and Nass clusters (see Fig. 0}, having nearly 
spherical geometries, conforms the jellium description. 
For instance, Nass shows a jellium-like behavior with 
s, p, d, ... shell structure. However, for systems such 
as Na25 and Naso, due to the disordered nature of the 
ground-state geometries, the degeneracy in the eigen- 
value spectra is lifted, thereby leading to a continuous 
spectra. 



2.0 



1.0 

(0 
0) 

I 2.0 

o 

? 1.0 

0) 
Q. 
CO 

■o 

0) 
N 



2.0 



1-0 

E 

^ 2.0 



Q 
. O 


. 1 u 






. 1 o 


1 o 






on 


25 










M\ 

■ 


50 












92 

J 


I 



1.0 
2.0 
1.0 

200 400 200 400 
Temperature (K) 

FIG. 5: The normalized specific-heat as a function of temper- 
ature. Co = (3A^ — 9/2)fcs is the zero-temperature classical 
limit of the rotational plus vibrational canonical specific-heat 



B. Thermodynamics 

The thermodynamic behavior is studied by analyzing 
several indicators. We have calculated the specific-heat, 
caloric curve, Lindemann criterion (Srms) as a function of 
temperature for each cluster. These are shown in Figs.jSl 
and [7| respectively. The melting temperatures, Tm, 
taken as the temperature corresponding to the peak in 
the specific-heat curve, are shown in Fig. (HI 

The examination of thermodynamic indicators as a 
function of cluster size reveals interesting trends. It may 
be seen that the clusters of sizes 8 and 10, do not show 
any recognizable peak in the heat capacities. This is 
also reflected in the caloric curves which increase contin- 
uously. While the peak is rather broad for n = 13-20, 
it progressively becomes narrower as the size increases, 
and at n = 92, a rather sharp peak with a width of the 
order of 30 K is observed. The caloric curve (Fig. 
and the Srms (Fig- Ej) for larger clusters, viz. n = 40, 
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FIG. 6: The averaged potential energy (caloric curve) as a 
function of temperature. 



55, and 92, clearly show distinct solidlike, liquidlike and 
transition regions. Interestingly, the melting tempera- 
tures depicted in Fig. [S] show irregular pattern with the 
maximum variation of about 60 K. High melting tem- 
peratures are observed in two clusters: Na4o and Nass, 
which are very symmetric; Na4o exhibiting electronic clo- 
sure and Na55 representing a geometrically closed system. 
We plot 5E in Fig.El which is an indicator of anharmonic- 
ity in the system. It is clear from the figure that for the 
clusters with n > 40, it is possible to distinguish a tem- 
perature region (<200 K) showing harmonic behavior. 
In contrast to this, the smaller clusters change from har- 
monic to anharmonic behavior nearly continuously. The 
most remarkable observation concerns the trends in the 
specific-heat curve and other indicators for sizes n — 40, 
50 and 55. In spite of a well-defined peak in the specific- 
heat curves for Na4o and Nass, the Naso cluster shows a 
rather broad structure, similar to that seen in smaller 
clusters. This size-sensitive behavior is discussed further 
below. 

Thus, all the indicators, like the specific-heat, the 
caloric curve, 5rms and 5E clearly show that small clus- 
ters (n = 8, 10) do not undergo any melting-like tran- 
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FIG. 7: The root-mean-square bond length fluctuation 
(<5rms) as a function of temperature. 
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FIG. 8: The melting temperature as a function of size. The 
melting temperature in bulk (Tm[Bulk]) is 370 K is shown 
with the horizontal line. 
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monic limit plotted as a function of temperature. The error 
bar is the standard error at each temperature. 



sition. The examination of their ionic motion indicates 
that over the entire range of the temperatures the motion 
is dominated by isomer hopping. This is in agreement 
with the observation by Manninen et al. However, 
these results are in contrast with the SMA and tight bind- 
ing calculations of Calvo et al. Q , and the Car-Parrinello 
orbital free calculations of Aguado et al. Recall that 
these calculations, though are not in agreement with each 
other, show distinct peaks (broader in case of SMA) in 
the heat capacities of these clusters. Our simulations fur- 
ther show the smaller clusters to be dissociated at about 
750 K. Our results for Naia and Na2o are also in dis- 
agreement with those by Calvo et al. 8] They calculated 
the canonical heat capacity of Naia, using an icosahe- 
dron for SMA calculations and a pentagonal structure 
with Ci symmetry for TB calculations as the ground- 
state geometries. While the heat capacity with the SMA 
potential exhibited a single prominent peak, that of TB 
calculation showed a premelting feature. They attributed 
this difference to the difference in the ground-state ge- 
ometries. They further found such premelting feature in 
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FIG. 10: The difference of average potential energy SEpot 
between solid-state and liquid-state as a function of size. 



the heat capacity of Na2o cluster, for which they used 
a capped double icosahedron as the lowest energy struc- 
ture. However, our calculations show that this structure 
is about 0.1 eV higher than the ground-state structure 
obtained by us for Na2o. Thus, the differences in the 
specific-heat may be also caused by the differences in 
the ground-state geometries. The Srms for all the clus- 
ters (Fig. [Tj) clearly shows that for smaller systems (n < 
40), it increases almost continuously, whereas for larger 
ones, a sharp rise is seen that corresponds to the peak 
in the specific-heat. The behavior of SE is, as noted 
earlier, consistent with this. We have also calculated a 
quantity SEpot defined as the difference of average po- 
tential energy of the melted cluster with respect to the 
ground-state structure at T = K. Schmidt et al.^E^ 
have inferred from the experimental caloric curve that 
the melting temperature is strongly influenced by such 
an energy contribution. They showed further that Tm 
follows closely the variation in energy difference between 
solid and liquid as a function of the cluster size. To ex- 
amine this feature, we plot the potential energy differ- 
ence between solidlike-state and liquidlike-state shown 
in Fig.^l For this purpose, we have taken the high tem- 
perature as 300 K for n = 13, 15, 20, 25, and 50, 350 K 
for n = 40, 55, which are about 50 K to 70 K higher 
than the melting temperature of the respective clusters. 
It can be immediately seen that the variation of SEpot 
follows that of the melting temperature (Fig. |SJ| , clearly 
indicating that the melting transition is mainly driven by 
energy contribution. 



C. Na4o • Naso • Nass 

Now, we turn our discussion to the most remarkable 
observation concerning the trends in the specific-heat in 
going from n = 40 to n = 55. As mentioned earlier, while 
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Na4o and Nass show well-defined peaks in the specific- 
heat curve (peak for Nass being much sharper), the peak 
for Naso is rather flat and is almost similar to that of the 
smaller clusters (say, n = 13-20). We note that Naso, 
being larger than Na4o, is expected to show slightly bet- 
ter melting transition. Interestingly, what is seen is ex- 
actly the opposite. It may be noted that such peculiar 
size-sensitivity has been observed experimentally in two 
systems, namely clusters of gallium (Ga+, n = 30-50 and 
55) and clusters of aluminum (A1+, n = 49-62) 0. 
For instance, in case of Ga clusters the heat capacity for 
n = 30 shows a flat curve without a peak. For n = 31 
it shows a remarkably sharp peak. Interestingly, addi- 
tion of one more atom {i.e. n — 32) diminishes this peak 
making the heat capacity nearly flat. We believe this 
behavior to be generic as it has not only been observed 
experimentally in case of gallium clusters but also for 
aluminum cluster and theoretically for sodium clusters 
in the present simulations. 

We note certain peculiar characteristics of the thermo- 
dynamic properties of Na4o, Naso, and Nass. We find the 
melting temperature of Naso to be about 60 K lower that 
that of Na4o and Nass . The Srms for Naso exhibits grad- 
ual increase in the temperature range of 100 K to 300 K. 
Further, the energy deviation SE for Naso, as seen in 
Fig. |51 starts to increase continuously at about 100 K 
to up to about 400 K, indicating a continuous change 
from harmonic behavior to anharmonic one. This behav- 
ior of Naso is in contrast with that of Na4o and Nass, 
where the change is seen in a smaller temperature width 
(of about 30-40 K) around the melting temperature. In 
order to bring out the origin of this phenomena, we exam- 
ine the nature of the ground-state geometries for these 
three clusters. Na4o and Nass are very symmetric struc- 
tures having almost five-fold symmetry. The values of 
the shape deformation parameter, edef, shown in Fig. O 
clearly indicate that they are nearly spherical structures. 
Further, the eigenvalue spectra of the ground-state ge- 
ometries of Na4o and Nass clusters also show this symme- 
try, conforming the jellium model. However, the eigen- 
value spectrum of Naso is very different from those of 
Na4o and Nass, in the sense that there are levels in the 
energy gaps leading to a more uniform spectrum. In this 
sense, Na4o and Nass are ordered, i.e. more symmetric, 
and Naso is amorphous. 

In Fig. ^2 we show the distances from the center of 
mass of all the atoms in the ground-state geometries of 
Naso and Nass ■ Clearly, the ordered geometric shell 
structure of Nass is destroyed when five atoms are re- 
moved. We believe that the nature of ground-state ge- 
ometry of the cluster has a significant effect on its melt- 
ing characteristics. An ordered or symmetric cluster, like 
Na4o and Nass , is expected to give rise to a well-defined 
peak in the heat capacity, while an amorphous and dis- 
ordered clusters, like Na2s and Naso, may lead to a con- 
tinuous melting transition. 
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FIG. 11: The distance from the center of mass for the 
ground-state geometries of Naso and Nass. 



IV. SUMMARY 

We have investigated the thermodynamics properties 
of small sodium clusters, Na„, in the size range of 8 to 55, 
using ab initio molecular dynamics with simulation time 
of 1.3-1.8 ns per cluster. We have analyzed several ther- 
modynamic indicators such as the specific-heat, caloric 
curve, Lindemann criterion, and the deviation energy SE 
to understand the melting characteristics in these clus- 
ters. We observe irregular variation in the melting tem- 
peratures as a function of size, which has also been seen 
in experiments by Haberland and coworkers for larger 
clusters. The reduction of about 30 % than the bulk 
value in melting temperature of sodium clusters, seen 
in experiment is also observed here. Further, we find 
a strong correlation between the ground-state geometry 
and the finite temperature characteristics of the sodium 
clusters. If a cluster has an ordered geometry, it is likely 
to show a relatively sharp melting transition. However, 
a cluster having a disordered geometry is expected to ex- 
hibit a broad peak in the specific-heat curve, indicating 
a poorly-defined melting transition. The size-sensitivity 
in the melting transition, seen in experiments by Breaux 
et al. 00, is observed for the case of sodium clusters in 
the size range of 40 to 50. The calculation of potential 
energy difference between solidlike state and liquidlike 
state reveals that the melting transition is mainly driven 
by the energy contribution and the entropy has a minor 
role in melting phenomena. 
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